#pragma once
#ifndef POKLINGTON

#define POKLINGTON
#include "../Lab2_Primality_test/solovay_strassen_test.h"
#include "../Lab1_Large_numbers/__uintX2.h"

#include <math.h>
#include <vector>
#include <assert.h>
#include <new>

typedef unsigned __int16 small_number;
typedef std::vector<uint128> small_list;

// Poklington primality test
class poklington
{
private:
	small_number* small_prime;
	small_number small_count;
	//Sieve of Eratosthenes
	void sieve(const small_number max)
	{
		if(max < 2) throw "Not enough prime numbers";
		bool* prime=NULL;
		try
		{
			prime=new bool[max+1];
			try
			{
				// 1 is prime, 2 is prime
				prime[0]=prime[1]=false;
				// all the rest consider composite
				for(small_number i=2; i<=max; ++i) prime[i]=true;


				for(small_number i=2; i<=ceil(sqrtl(max)); ++i) 
				{
					if(i > max) throw "Index out of range";
					if(prime[i])
					{
						for(small_number j=i; j<= max/i; ++j) prime[i*j]=false;
					}
				}

				for(small_number i=2; i<=max; ++i) 
					if(prime[i]) ++small_count;

				small_prime=new small_number[small_count];
		
				for(small_number i=0, j=0; i<=max; ++i) 
					if(prime[i]) small_prime[j++]=i;
			}
			catch(char* e)
			{
				printf(e);
			}
			delete[] prime;
		}
		catch (std::bad_alloc)
		{
			if(NULL != prime) 
				delete[] prime;
		}
	}
	static bool is_prime(const small_list factors,
			   //const usize_t fact_length,
			   const uint128 R,
			   const unsigned __int8 t, uint128 &g=uint128(0)) //const
	{
		uint128 n=R;
		for(usize_t i = 0; i < factors.size(); ++i)
		{
			n *= factors[i];
		}
		++n;
		/*  N = R*F + 1  */
		//if((err>=1)||(err<=0)) throw "Error possibility should be between 0 and 1";
		//const unsigned __int8 t=(int)ceil(log(1/err)/log(2.0));
		//if((n&1)==0) return false;
		for(unsigned __int8 iter=0;iter<t;++iter)
		{
			uint128 a=uint128::random()%(n-1)+1;
			if(uint128::pow_m(a,n-1,n)!=1) 
				return false; /* not satysfied Ferma's theorem result,
							     definetely composite  */
			
			bool composite=false;
			for(unsigned __int8 j=0;j<factors.size();++j)
			{
				bool repeat=false;
				for(unsigned __int8 k=0;k<j;++k)
				{
					if(factors[j]==factors[k]) //already been?
					{
						repeat=true;
						break;
					}
				}
				if(repeat) continue;

				if(uint128::pow_m(a,(n-1)/factors[j],n)==1) // th.Selfridge
				{
					composite=true;
					break;
				}
			}
			if(!composite) 
			{
				g=a;
				return true; //definetely prime
			}
		}
		return false; //probably composite
	}


public:
	poklington(const unsigned __int16 MAX=500):small_count(0){sieve(MAX);}
	~poklington(){delete[] small_prime;}

	unsigned int prime_attempts;
//#define SCREEN_LOG

	uint128 get_prime(const usize_t N_SIZE, uint128 &gen=uint128(0))
	{
		return get_prime(1, N_SIZE, gen);
	}
	// find and return N-bit prime number p and (optionally) generator of set Z_p*
	uint128 get_prime(const uint128 &start_prime, const usize_t N_SIZE, uint128 &gen=uint128(0))
	{
		if(N_SIZE > sizeof(uint128)*8) 
			throw "The size cannot be greater than size of maximum available value";
		
		prime_attempts=0;
		
		while(true)
		{
			small_list factors_list;
			uint128 F(start_prime);

			if(start_prime != uint128::one())
				factors_list.push_back(start_prime);
			
			while(F.size()<=N_SIZE/2) // F should be greater than R (th. Poklington)
			{
				small_number factor=small_prime[rand()%small_count];
				factors_list.push_back(factor);
				F*=factor; // F contains known small factors
			}
			
			if(N_SIZE <= F.size()) continue;
			
			usize_t R_SIZE=N_SIZE-F.size();
			if(R_SIZE < 2) continue; 

			uint128 R=((uint128::random() & 
				       ((uint128(1)<<(R_SIZE-2))-1))//set R_SIZE-2 random bits
				              <<1)              //set last zero-bit to 0, cause R should be even
					 |(uint128(1)<<(R_SIZE-1)); //set first R_SIZE-1-th bit to 1 cause R should be length of R_SIZE
            
			uint128 N = R*F + 1;
			if(N.size()!=N_SIZE) continue;

#ifdef SCREEN_LOG
			std::cout<<"probing "<<N<<" with size of "<<N.size()<<"\n";
#endif /* SCREEN_LOG */
			/*
			small_number fact_length=factors_list.size();
			small_number* factors=new small_number[fact_length];
			for(small_number i=0;i<fact_length;++i)
			{
				factors[i]=factors_list.front(); //generating array from queue of prime small factors
				factors_list.pop();
			}
			*/
			if(is_prime(factors_list,R,2,gen)) 
			{
				assert(probably_prime(N)); 
				return N;
			}

			//delete[] factors;
			if(probably_prime(N)) 
			{
				++prime_attempts;
#ifdef SCREEN_LOG
				std::cout<<"\tprobably composite\n";
#endif /* SCREEN_LOG */
			}
			else
			{
#ifdef SCREEN_LOG
				std::cout<<"\tdefinitely composite\n";
#endif /* SCREEN_LOG */
			}
		}
		return 0;
	}
};
#endif /* POKLINGTON */